home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / LIBIP / FFT.C < prev    next >
C/C++ Source or Header  |  1999-09-11  |  3KB  |  86 lines
  1. /* 
  2.  * fft.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. #define PI    3.14159265358979
  10.  
  11. /*
  12.  * fft()
  13.  *   DESCRIPTION:
  14.  *       Uses time decomposition with input bit reversal. The Cooley/Tookey
  15.  *       Fortran scheme for doing recursive odd/even decimation without really
  16.  *       doing bit reversal is used. The computation is done in place, so the
  17.  *       output replaces the input.  The contents of the arrays are changed
  18.  *     from the input data to the FFT coefficients.
  19.  *   ARGUMENTS:
  20.  *       npoints(int) the number of points in the FFT. Must be a power of 2.
  21.  *       real, imag(float *) pointers to arrays of floats for input and output.
  22.  *       inv(int)       inv = 1 for inverse transform
  23.  *                      inv = -1 for forward transform
  24.  *   RETURN VALUE:
  25.  *     none
  26.  */
  27. void
  28. fft (npoints, real, imag, inv)
  29.      int npoints, inv;
  30.      float *real, *imag;
  31. {
  32.   register int i, index, swapindex, j, k;
  33.   register float tr, ti, angle, wr, wi;
  34.   double sin (), cos ();
  35.  
  36.   /* SWAP THE INPUT ELEMENTS FOR THE DECIMATION IN TIME ALGORITHM. */
  37.   for (index = 1, swapindex = 0; index < npoints; index++) {
  38.     k = npoints;
  39.     do
  40.       k /= 2;
  41.     while ((swapindex + k) >= npoints);
  42.     swapindex = (swapindex % k) + k;
  43.     if (swapindex <= index)
  44.       continue;
  45.     tr = real[index];
  46.     real[index] = real[swapindex];
  47.     real[swapindex] = tr;
  48.     ti = imag[index];
  49.     imag[index] = imag[swapindex];
  50.     imag[swapindex] = ti;
  51.   }
  52.  
  53.   /*
  54.    * DO THE BUTTERFLY COMPUTATIONS.
  55.    *      stage index:    k = 1, 2, 4, 8 . . .
  56.    *                      For npoints = 8, for example, there will
  57.    *                      be three stages of butterfly computations.
  58.    *      b'fly indices:  i and j will be separated by a distance
  59.    *                      dependent on the current stage.
  60.    *                      k is used as the separation constant.
  61.    */
  62.  
  63.   for (k = 1; k < npoints; k *= 2)
  64.     for (index = 0; index < k; index++) {
  65.       angle = (float) PI *((float) index * inv) / ((float) k);
  66.       wr = (float) cos (angle);
  67.       wi = (float) sin (angle);
  68.       for (i = index; i < npoints; i += 2 * k) {
  69.         j = i + k;
  70.         tr = (wr * (real[j])) - (wi * (imag[j]));
  71.         ti = (wr * (imag[j])) + (wi * (real[j]));
  72.         real[j] = real[i] - tr;
  73.         imag[j] = imag[i] - ti;
  74.         real[i] += tr;
  75.         imag[i] += ti;
  76.       }
  77.     }
  78.  
  79.   /* for inverse transform, scale output by 1/N */
  80.   if (inv == 1)
  81.     for (i = 0; i < npoints; i++) {
  82.       real[i] = real[i] / npoints;
  83.       imag[i] = imag[i] / npoints;
  84.     }
  85. }
  86.